# Remove columns
air_pollution_df_cleaned <- subset(air_pollution_df, select = -c(Indoor.air.pollution..deaths.per.100.000., Outdoor.particulate.matter..deaths.per.100.000., Outdoor.ozone.pollution..deaths.per.100.000.))

# Set datatypes
air_pollution_df_cleaned$Entity = factor(air_pollution_df_cleaned$Entity)
air_pollution_df_cleaned$Code = factor(air_pollution_df_cleaned$Code)
air_pollution_df_cleaned$Year = factor(air_pollution_df_cleaned$Year)

# Rename columns
air_pollution_df_cleaned <- rename(air_pollution_df_cleaned, Country = Entity, ISO.alpha3.code = Code, Deaths.Air.Pollution.per.100k = Air.pollution..total...deaths.per.100.000.)

str(air_pollution_df_cleaned)

# Remove columns
gdp_df_cleaned <- subset(gdp_df, select = -c(Indicator.Name, Indicator.Code))

# Pivot wide to long dataset (data melt)
gdp_df_cleaned <- gdp_df_cleaned %>%
  pivot_longer(
    cols = starts_with("X"), 
    names_to = "Year", 
    values_to = "GDP.USD", 
    values_drop_na = TRUE
  )

# Remove characters from column
gdp_df_cleaned$Year<-gsub("X","",as.character(gdp_df_cleaned$Year))

# Set datatypes
gdp_df_cleaned$Country.Name = factor(gdp_df_cleaned$Country.Name)
gdp_df_cleaned$Country.Code = factor(gdp_df_cleaned$Country.Code)
gdp_df_cleaned$Year = factor(gdp_df_cleaned$Year)

# Rename columns
gdp_df_cleaned <- rename(gdp_df_cleaned, Country = Country.Name, ISO.alpha3.code = Country.Code)

str(gdp_df_cleaned)

# Remove columns
population_region_df_cleaned <- subset(population_region_df, select = -c(Notes, Type, Parent.code))

# Pivot wide to long dataset (data melt)
population_region_df_cleaned <- population_region_df_cleaned %>%
  pivot_longer(
    cols = starts_with("X"), 
    names_to = "Year", 
    values_to = "Population.thousands", 
    values_drop_na = TRUE
  )

# Remove characters from column
population_region_df_cleaned$Year<-gsub("X","",as.character(population_region_df_cleaned$Year))
population_region_df_cleaned <- as.data.frame(apply(population_region_df_cleaned, 2, function(x) gsub("\\s+", "", x))) 

# Set datatypes
population_region_df_cleaned$Country = factor(population_region_df_cleaned$Country)
population_region_df_cleaned$SDGRegion = factor(population_region_df_cleaned$SDGRegion)
population_region_df_cleaned$Country.code = factor(population_region_df_cleaned$Country.code)
population_region_df_cleaned$SubRegion = factor(population_region_df_cleaned$SubRegion)
population_region_df_cleaned$Year = factor(population_region_df_cleaned$Year)
population_region_df_cleaned$Population.thousands = as.numeric(population_region_df_cleaned$Population.thousands)

# Rename columns
population_region_df_cleaned <- rename(population_region_df_cleaned, M49.code = Country.code)

str(population_region_df_cleaned)

# Set datatypes
country_codes_df$M49.code = factor(country_codes_df$M49.code)
country_codes_df$Country.or.Area = factor(country_codes_df$Country.or.Area)
country_codes_df$ISO.alpha3.code = factor(country_codes_df$ISO.alpha3.code)
country_codes_df$ISO.alpha2.code = factor(country_codes_df$ISO.alpha2.code)

str(country_codes_df)

# Set datatypes
world$iso_a2 = factor(world$iso_a2)

# Remove columns
world_df_cleaned <- subset(world, select = c(iso_a2, geom))

str(world_df_cleaned)

# clean Air Quality Dataset for merge
airquality <- read.csv("/Users/b/Desktop/DataScience/DataScience/GitHub_Repo/Midterm Project/data/airquality.csv")
str(airquality)

airquality$ISO.alpha3.code<-airquality$LOCATION
airquality$Pollution.Exposure<-airquality$SUBJECT

airquality <- subset(airquality, select = c(ISO.alpha3.code, Year, AirQ,Pollution.Exposure))
str(airquality)

# Join datasets
final_df <- merge(x = air_pollution_df_cleaned, y = country_codes_df, by = 'ISO.alpha3.code')
final_df <- merge(x = final_df, y = gdp_df_cleaned, by = c('ISO.alpha3.code', 'Year'))
final_df <- merge(x = final_df, y = population_region_df_cleaned, by = c('M49.code', 'Year'))
final_df <- merge(x = final_df, y = world_df_cleaned, by.x = 'ISO.alpha2.code', by.y = 'iso_a2', all.x = TRUE) #left join
final_df <- merge(x= final_df, y= airquality, by=c("Year","ISO.alpha3.code")) 


# Remove columns
final_df <- subset(final_df, select = -c(Country.or.Area, Country.y, Country))

# Calculate GDP per capita 
final_df$gdp.per.capita <- final_df$GDP.USD / final_df$Population.thousands

str(final_df)

#Deaths
death_hist <- ggplot() + 
      geom_histogram(data = final_df, aes(x = Deaths.Air.Pollution.per.100k), bins=30, color="darkblue", fill="steelblue", binwidth = 5)  +
              labs(title = "Distribution of Deaths per 100,000 from Air Pollution",
                          subtitle = "X Countries, 1990 - 2017", 
                          caption = "Data Source: Kaggle",
                          y = "Count",
                          x = "Deaths per 100,000") +
        theme_minimal() +
        theme(axis.title = element_text(size = 8, face = "bold"),
              panel.grid.major.x = element_blank(),
              panel.grid.minor = element_blank(),
              panel.background = element_blank(),
              panel.grid.major.y = element_blank(),
              axis.line.x = element_line(color = "black"),
              axis.ticks = element_line(color = "black"),
              axis.text = element_text(size = 10),
              legend.position = "none",
              plot.subtitle = element_text(size = 8),
              plot.title = element_text(size = 10, margin = margin(b = 10)))

death_hist <- death_hist + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                           plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                           plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))

death_hist

##Population
pop_hist <- ggplot() + geom_histogram(data = final_df, aes(x = Population.thousands), bins=30, color="darkgreen", fill="green4")  +
              labs(title = "Population Distribution",
                          subtitle = "X Countries, 1990 - 2017", 
                          caption = "Data Source: United Nations",
                          y = "Count",
                          x = "Population") +
              theme_minimal() +
              theme(axis.title = element_text(size = 8, face = "bold"),
                    panel.grid.major.x = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    panel.grid.major.y = element_blank(),
                    axis.line.x = element_line(color = "black"),
                    axis.ticks = element_line(color = "black"),
                    axis.text = element_text(size = 10),
                    legend.position = "none",
                    plot.subtitle = element_text(size = 8),
                    plot.title = element_text(size = 10, margin = margin(b = 10)))

pop_hist <- pop_hist + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                           plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                           plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))

pop_hist

##GDP
gdp_hist <- ggplot() + geom_histogram(data = final_df, aes(x = gdp.per.capita), bins=30, color="orangered4", fill="orangered3")  +
              labs(title = "Distribution of GDP per Capita",
                          subtitle = "X Countries, 1990 - 2017", 
                          caption = "Data Source: United Nations",
                          y = "Count",
                          x = "GDP per Capita") +
              theme_minimal() +
              theme(axis.title = element_text(size = 8, face = "bold"),
                    panel.grid.major.x = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    panel.grid.major.y = element_blank(),
                    axis.line.x = element_line(color = "black"),
                    axis.ticks = element_line(color = "black"),
                    axis.text = element_text(size = 10),
                    legend.position = "none",
                    plot.subtitle = element_text(size = 8),
                    plot.title = element_text(size = 10, margin = margin(b = 10)))

gdp_hist <- gdp_hist + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                           plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                           plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))

gdp_hist

death_boxplot <- ggplot(final_df, aes(x=SDGRegion, y = Deaths.Air.Pollution.per.100k)) + 
                  geom_boxplot(fill="skyblue")
death_boxplot <- death_boxplot + labs(title="BoxPlot of Deaths per 100,000 from Air Pollution, by SDG Region", x = "SDG Region", y ="Deaths per 100k from Air Polution") + 
  theme_classic() + theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 5))
death_boxplot

sub_boxplot <- ggplot(final_df, aes(x=SubRegion, y = Deaths.Air.Pollution.per.100k)) + 
                  geom_boxplot(fill="azure3")
sub_boxplot <- sub_boxplot + labs(title="BoxPlot of Deaths per 100k from Air Pollution, by Sub Regions", x = "Sub Regions", y ="Deaths per 100k from Air Polution") + 
  theme_classic() + theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 4))
sub_boxplot

gdp_boxplot <- ggplot(final_df, aes(x=SubRegion, y = gdp.per.capita)) + 
                  geom_boxplot(fill="aquamarine3")
gdp_boxplot <- gdp_boxplot + labs(title="BoxPlot of GDP per Capita, by Sub Regions", x = "Sub Regions", y ="GDP per Capita (USD per person)") + 
  theme_classic() + theme(axis.text.x = element_text(angle = 30, hjust = 1, size = 4))
gdp_boxplot

# Convert to sf object so we can plot it
final_df_sf = st_as_sf(final_df)

# Choose relevant columns
final_df_sf_regions <- subset(final_df_sf, select = c(SDGRegion, SubRegion))
final_df_sf_intensities <- subset(final_df_sf, select = c(gdp.per.capita, Deaths.Air.Pollution.per.100k, Population.thousands))

plot(final_df_sf_regions)

plot(final_df_sf_intensities)

g1 = ggplot() + geom_sf(data = final_df_sf, aes(fill = Deaths.Air.Pollution.per.100k), color = 'black') + scale_fill_distiller(palette = "YlGnBu") + coord_sf(crs = st_crs(4283)) + theme_minimal()

g1

g2 = ggplot() + geom_sf(data = final_df_sf[final_df_sf$SDGRegion == 'EASTERNANDSOUTH-EASTERNASIA' & final_df_sf$Year == 2017, ], aes(fill = Deaths.Air.Pollution.per.100k), color = 'black') + scale_fill_distiller(palette = "YlGnBu") + coord_sf(crs = st_crs(4283)) + theme_minimal()

g2

str(final_df)

qqnorm(final_df$Population.thousands)

qqnorm(final_df$AirQ)

cor(final_df$Population.thousands, final_df$AirQ, method = c("pearson"))

cor.test(final_df$Population.thousands, final_df$AirQ, method=c("pearson"))

loadPkg("ggpubr")
ggscatter(final_df, x = "Population.thousands", y = "AirQ", 
          add = "reg.line", conf.int = TRUE, 
          cor.coef = TRUE, cor.method = "pearson",
          xlab = "Population in thousands", ylab = "Pollution Exposure")

# Plot and Regression line
plot(final_df$Population.thousands, final_df$AirQ, pch = 19, col = "navyblue")
abline(lm(final_df$AirQ ~ final_df$Population.thousands), col = "red", lwd = 3)
text(paste("Correlation:", round(cor(final_df$Population.thousands, final_df$AirQ), 2)), x = 25, y = 95)

# Pearson correlation
pop_poll_cor<- cor(select(final_df, AirQ, Population.thousands, gdp.per.capita,))
pop_poll_cor

#correlation matrix
xkabledply(pop_poll_cor)

#plot
loadPkg("corrplot")
corrplot(pop_poll_cor)

library("data.table")
lowairq <- subset(final_df, AirQ >51 & Year == 2017)
str(lowairq)

lowairq_top <- lowairq[order(lowairq$AirQ, decreasing = TRUE), ]  # Top N highest values by group
lowairq_top <- data.table(lowairq_top, key = "Country.x")
lowairq_top <- lowairq_top[ , head(.SD, 50)]
str(lowairq_top)

air1 <- ggplot(data = lowairq_top,
             mapping = aes(x = AirQ, y = reorder(Country.x , AirQ), color = SDGRegion)) +
      geom_segment(data = lowairq_top,
               mapping = aes(x = 0, xend = AirQ, y = reorder(Country.x, AirQ), yend = Country.x)) +
      geom_point() +
          scale_x_continuous(limits = c(0, 150), breaks = seq(0,150,25)) +
          scale_color_manual(values = c("darkorange3", "firebrick1","lightpink1","green3","dodgerblue4","gold1")) 

air1 <- air1 + labs(title = "Top 50 Highest Pollution Exposure, by Country and Region",
                            subtitle = "Average Exposure to PM2.5 fine particles, 2017", 
                            caption = "Data Source: Air quality and health: Exposure to PM2.5 fine particles (OECD)",
                            x = "Pollution Exposure",
                            y = "") +
         theme_minimal() +
        theme(axis.title = element_text(size = 8, face = "bold"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.ticks = element_line(color = "black"),
        axis.text = element_text(size = 10),
        #legend.position = "none",
        plot.subtitle = element_text(size = 8),
        plot.title = element_text(size = 10, margin = margin(b = 10)))

air1 <- air1 + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                           plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                           plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))
air1

highairq <-subset(final_df, AirQ <51 & Year == 2017)
str(highairq)

highairq_top <- highairq[order(highairq$AirQ, decreasing = FALSE), ]  # Top N lowest values by group
highairq_top <- data.table(highairq_top, key = "Country.x")
highairq_top <- highairq_top[ , head(.SD, 50)]
str(highairq_top)

air2 <- ggplot(data = highairq_top,
             mapping = aes(x = AirQ, y = reorder(Country.x , AirQ), color = SDGRegion)) +
      geom_segment(data = highairq_top,
               mapping = aes(x = 0, xend = AirQ, y = reorder(Country.x, AirQ), yend = Country.x)) +
      geom_point() +
          scale_x_continuous(limits = c(0, 150), breaks = seq(0,150,25)) +
          scale_color_manual(values = c("darkorange3", "firebrick1","lightpink1","green3","dodgerblue4","gold1", "mediumorchid1", "azure4","sienna3")) 

air2 <- air2 + labs(title = "Top 50 Lowest Pollution Exposure, by Country and Region",
                            subtitle = "Average Exposure to PM2.5 fine particles, 2017", 
                            caption = "Data Source: Air quality and health: Exposure to PM2.5 fine particles (OECD)",
                            x = "Pollution Exposure",
                            y = "") +
         theme_minimal() +
        theme(axis.title = element_text(size = 8, face = "bold"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.ticks = element_line(color = "black"),
        axis.text = element_text(size = 10),
        #legend.position = "none",
        plot.subtitle = element_text(size = 8),
        plot.title = element_text(size = 10, margin = margin(b = 10)))

air2 <- air2 + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                           plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                           plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))
air2

library(scales)
region <- ggplot(final_df[which(final_df$Year ==2017),], aes(y= reorder(SDGRegion,AirQ), x= AirQ, fill = "dodgerblue4")) +
  geom_bar(stat ="identity", size = 2.5, width=0.5, position = position_dodge(width=0.5)) +
  scale_fill_manual(values=c("dodgerblue4")) +
  scale_x_continuous(label=comma, limits=c(0,150), breaks = seq(0,150,25))

region <- region + labs(title = "Pollution Exposure, by Region",
                                          subtitle = "Average Exposure to PM2.5 fine particles, 2017", 
                                          caption = "Air quality and health: Exposure to PM2.5 fine particles (OECD)",
                                          y = "",
                                          x = "Pollution Exposure") +
  theme_minimal() +
  theme(axis.title = element_text(size = 8, face = "bold"),
        panel.grid.major.x = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.grid.major.y = element_blank(),
        axis.line.x = element_line(color = "black"),
        axis.ticks = element_line(color = "black"),
        axis.text = element_text(size = 10),
        legend.position = "none",
        plot.subtitle = element_text(size = 8),
        plot.title = element_text(size = 10, margin = margin(b = 10))) +
  scale_color_manual(values = c("darkblue"))


region <- region + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                                           plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                                           plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))
region

library(scales)
#Aggregate data by total deaths by year
deathsbyyear <- group_by(.data = final_df,Year)
totdeath <- summarize(.data = deathsbyyear, tot_deaths = sum(Deaths.Air.Pollution.per.100k, na.rm = TRUE))
str(totdeath)

deaths_line <- ggplot() +
               geom_line(data = totdeath, mapping = aes(x = Year, y = tot_deaths, group = 1), size = 1.2) +
              geom_point(data = totdeath, mapping = aes(x = Year, y = tot_deaths, group = 1), size = 1.2) +
              scale_color_manual(values = c("darkmagenta")) +
              scale_y_continuous(label = comma, limits = c(0, 40000), breaks = seq(0,40000,10000)) 


deaths_line <- deaths_line + labs(title = "Deaths per 1000,000 by Air Pollution",
                                            subtitle = "1990 - 2017", 
                                            caption = "?",
                                            y = "air Pollution Induced Deaths",
                                            x = "") +
              theme_minimal() +
              theme(axis.title = element_text(size = 8, face = "bold"),
                    panel.grid.major.x = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    panel.grid.major.y = element_blank(),
                    axis.line.x = element_line(color = "black"),
                    axis.ticks = element_line(color = "black"),
                    axis.text = element_text(size = 10),
                    legend.position = "none",
                    plot.subtitle = element_text(size = 8),
                    plot.title = element_text(size = 10, margin = margin(b = 10)))
            
            deaths_line <- deaths_line + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                                         plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                                         plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))
            deaths_line

library(scales)
#Aggregate data by total deaths by year
deathsbyyear_reg <- group_by(.data = final_df,Year, SDGRegion)
totdeath_reg <- summarize(.data = deathsbyyear_reg, total = sum(Deaths.Air.Pollution.per.100k, na.rm = TRUE))
str(totdeath)

deaths_line <- ggplot() +
               geom_line(data = totdeath_reg, mapping = aes(x = Year, y = total, group = SDGRegion, color = SDGRegion), size = 1.2) +
              geom_point(data = totdeath_reg, mapping = aes(x = Year, y = total, group = SDGRegion, color = SDGRegion), size = 1.2) +
              scale_y_continuous(label = comma, limits = c(0, 20000), breaks = seq(0,20000,2000)) 


deaths_line <- deaths_line + labs(title = "Deaths per 1000,000 by Air Pollution, by Region",
                                            subtitle = "1990 - 2017", 
                                            caption = "?",
                                            y = "air Pollution Induced Deaths",
                                            x = "") +
              theme_minimal() +
              theme(axis.title = element_text(size = 8, face = "bold"),
                    panel.grid.major.x = element_blank(),
                    panel.grid.minor = element_blank(),
                    panel.background = element_blank(),
                    panel.grid.major.y = element_blank(),
                    axis.line.x = element_line(color = "black"),
                    axis.ticks = element_line(color = "black"),
                    axis.text = element_text(size = 10),
                    #legend.position = "none",
                    legend.text = element_text(size=5),
                    plot.subtitle = element_text(size = 8),
                    plot.title = element_text(size = 10, margin = margin(b = 10)))
            
            deaths_line <- deaths_line + theme(plot.title = element_text(color = "black", size = 12, face = "bold", hjust = 0),
                                         plot.subtitle = element_text(color = "black", size = 10, hjust = 0 ),
                                         plot.caption = element_text(color = "black", size =8, face = "italic", hjust =0))
            deaths_line

str(air_pollution_df)

# select what we need
air_pollution_dfneed <- subset(air_pollution_df, select = -c(Code))

# make Entity and year factor
air_pollution_dfneed$Entity = factor(air_pollution_dfneed$Entity)
air_pollution_dfneed$Year = factor(air_pollution_dfneed$Year)

# sum up the death of each year
sum_total <- aggregate(air_pollution_dfneed$Air.pollution..total...deaths.per.100.000. , by=list(Year=air_pollution_dfneed$Year), FUN=sum)

names(sum_total)[names(sum_total) == 'x'] <- 'air_pollution_total'

sum_indoor <- aggregate(air_pollution_dfneed$Indoor.air.pollution..deaths.per.100.000. , by=list(Year=air_pollution_dfneed$Year), FUN=sum)

names(sum_indoor)[names(sum_indoor) == 'x'] <- 'air_pollution_indoor'

sum_outdoorP <- aggregate(air_pollution_dfneed$Outdoor.particulate.matter..deaths.per.100.000. , by=list(Year=air_pollution_dfneed$Year), FUN=sum)

names(sum_outdoorP)[names(sum_outdoorP) == 'x'] <- 'air_pollution_outdoorP'

sum_outdoorO <- aggregate(air_pollution_dfneed$Outdoor.ozone.pollution..deaths.per.100.000. , by=list(Year=air_pollution_dfneed$Year), FUN=sum)

names(sum_outdoorO)[names(sum_outdoorO) == 'x'] <- 'air_pollution_outdoorO'

# Merge the new variables by year
total1 <- merge(sum_total,sum_indoor, by="Year")
total2 <- merge(sum_outdoorO,sum_outdoorP, by="Year")
total3 <- merge(total1,total2, by="Year")

# Make graph using year as X with total death
library(ggplot2)
# total
ggplot(total3, aes(x=Year, y=air_pollution_total)) +
    geom_bar(stat='identity', position='dodge')

# indoor
ggplot(total3, aes(x=Year, y=air_pollution_indoor)) +
    geom_bar(stat='identity', position='dodge')

# outdoor particulate
ggplot(total3, aes(x=Year, y=air_pollution_outdoorP)) +
    geom_bar(stat='identity', position='dodge')

# outdoor ozone 
ggplot(total3, aes(x=Year, y=air_pollution_outdoorO)) +
    geom_bar(stat='identity', position='dodge')

# The max and min years
special1.1 <- sum_total[which.max(sum_total$air_pollution_total),]
special1.2 <- sum_total[which.min(sum_total$air_pollution_total),]
special1 <- rbind(special1.1, special1.2)  # Max and min of total

special2.1 <- sum_indoor[which.max(sum_indoor$air_pollution_indoor),]
special2.2 <- sum_indoor[which.min(sum_indoor$air_pollution_indoor),]
special2 <- rbind(special2.1, special2.2)  # Max and min of indoor

special3.1 <- sum_outdoorP[which.max(sum_outdoorP$air_pollution_outdoorP),]
special3.2 <- sum_outdoorP[which.min(sum_outdoorP$air_pollution_outdoorP),]
special3 <- rbind(special3.1, special3.2)  # Max and min of ondoorP

special4.1 <- sum_outdoorO[which.max(sum_outdoorO$air_pollution_outdoorO),]
special4.2 <- sum_outdoorO[which.min(sum_outdoorO$air_pollution_outdoorO),]
special4 <- rbind(special4.1, special4.2)  # Max and min of ondoorO

library("ggplot2")

# Year of max and min of total pollution
ggplot(special1, aes(x=Year, y=air_pollution_total)) +
    geom_bar(stat='identity', position='dodge')

# Year of max and min of indoor pollution
ggplot(special2, aes(x=Year, y=air_pollution_indoor)) +
    geom_bar(stat='identity', position='dodge')

# Year of max and min of ourdoor particulate pollution
ggplot(special3, aes(x=Year, y=air_pollution_outdoorP)) +
    geom_bar(stat='identity', position='dodge')

# Year of max and min of outdoor ozone pollution
ggplot(special4, aes(x=Year, y=air_pollution_outdoorO)) +
    geom_bar(stat='identity', position='dodge')

fit1 <- lm(Deaths.Air.Pollution.per.100k ~ gdp.per.capita, data=final_df_sf)

ggplot(final_df_sf, aes(gdp.per.capita, Deaths.Air.Pollution.per.100k)) +
  geom_point() +
  geom_smooth(method='lm', se=FALSE) +
  stat_regline_equation(label.y = 400, aes(label = ..eq.label..)) +
  stat_regline_equation(label.y = 350, aes(label = ..rr.label..))

final_df_sf_most_recent_5_yrs <- subset(final_df_sf, Year == c(2017, 2016, 2015, 2014, 2013))

fit2 <- lm(Deaths.Air.Pollution.per.100k ~ gdp.per.capita, data=final_df_sf_most_recent_5_yrs)

summary(fit2)
print(vif(fit2))

# plot(fit1)

ggplot(final_df_sf_most_recent_5_yrs, aes(gdp.per.capita, Deaths.Air.Pollution.per.100k)) +
  geom_point() +
  geom_smooth(method='lm')

final_df_sf_asia <- subset(final_df_sf, SDGRegion == 'EASTERNANDSOUTH-EASTERNASIA')

fit3 <- lm(Deaths.Air.Pollution.per.100k ~ log(gdp.per.capita), data=final_df_sf_asia)

summary(fit3)
print(vif(fit3))

# plot(fit1)

ggplot(final_df_sf, aes(log(gdp.per.capita), Deaths.Air.Pollution.per.100k)) +
  geom_point() +
  geom_smooth(method='lm') + 
  facet_wrap(~SDGRegion, ncol = 4)